function plotrocks(filenumber)
    N=readdata(0,1);
%     S=readdata(0,2);
%     T=readdata(0,3);
%     multigridlevels=readdata(0,4);
%     bound=readdata(0,5);
time=readdata(0,6);
%     dtnom=readdata(0,7);
dx=readdata(0,8);
dz=readdata(0,9);
mx=readdata(0,10);
mz=readdata(0,11);
%     mcs=readdata(0,12);
%     mct=readdata(0,13);
%     mtype=readdata(0,14);
%     mTemp=readdata(0,15);
%     mstressxx=readdata(0,16);
%     mstressxz=readdata(0,17);
%     mstrainxx=readdata(0,18);
%     mstrainxz=readdata(0,19);
%     mfinitestrain=readdata(0,20);
%     vx=readdata(0,21);
%     vz=readdata(0,22);
%     P=readdata(0,23);
     Nsurfaces=readdata(0,24);
     surfaces=readdata(0,25);

mcs=zeros(1,N); mct=zeros(1,N);
xdivision=1; xplotcellsize=sum(dx)/xdivision;
x=0; for s=0:xdivision-1 ind=find(mx>=x & mx<x+xplotcellsize); mcs(ind)=s; x=x+xplotcellsize; end
zdivision=30; zplotcellsize=sum(dz)/zdivision;
z=0; for t=0:zdivision-1 ind=find(mz>=z & mz<z+zplotcellsize); mct(ind)=t; z=z+zplotcellsize; end
Ia=int32(find((-1).^mcs.*(-1).^mct>0));
Ib=int32(find((-1).^mcs.*(-1).^mct<=0));
clear('mcs','mct');

mtotalstrain=zeros(N,1); oldtime=time;

figure(2); clf; contstdhist=[]; conthist=[]; timehist=[];
figure(3); clf; hold on;
for count=filenumber
    
    N=readdata(count,1);
    S=readdata(count,2);
    T=readdata(count,3);
    %     multigridlevels=readdata(count,4);
    %     bound=readdata(count,5);
    time=readdata(count,6);
    %     dtnom=readdata(count,7);
    dx=readdata(count,8);
    dz=readdata(count,9);
    mx=readdata(count,10);
    mz=readdata(count,11);
    %     mcs=readdata(count,12);
    %     mct=readdata(count,13);
    mtype=readdata(count,14);
    %mTemp=readdata(count,15);
    %mstressxx=readdata(count,16);
    %mstressxz=readdata(count,17);
    %mstrainxx=readdata(count,18);
    %mstrainxz=readdata(count,19);
    %mfinitestrain=readdata(count,20);
    %     vx=readdata(count,21);
    %     vz=readdata(count,22);
    %P=readdata(count,23);
    Nsurfaces=readdata(count,24);
    surfaces=readdata(count,25);

    dt=time-oldtime; oldtime=time;
    x=zeros(1,S); z=zeros(1,T);
    for s=1:S-1 x(s+1)=x(s)+dx(s); end
    for t=1:T-1 z(t+1)=z(t)+dz(t); end
    figure(1); clf;

    avtemp=zeros(1,T);

    mTemp=readdata(count,15);
    mk=zeros(1,N);
    mk(find(mtype==1))=3; mk(find(mtype==2))=3; mk(find(mtype>2))=4.08*power((298./(mTemp(find(mtype>2))+273)),0.406);
    [gridk,xgrid,zgrid]=bilingrid(mk,mx,mz,dx,dz,1,1);
    [gridtemp,xgrid,zgrid]=bilingrid(mTemp,mx,mz,dx,dz,1,1);
    gridq=zeros(S,T);
    for s=1:S
        for t=2:T-1
            gridq(s,t)=gridk(s,t)*(gridtemp(s,t+1)-gridtemp(s,t-1))/(dz(t-1)+dz(t));
        end
    end
    [Zgrid,Xgrid] = meshgrid(zgrid,xgrid); qsurface=interp2(Xgrid',Zgrid',gridq',x,surfaces(:,3)');

    mrho=zeros(1,N)+1;
    ind=find(mtype==1); mrho(ind)=2800-3.28e-5*mTemp(ind)*2800; ind=find(mtype==2); mrho(ind)=3000-3.28e-5*mTemp(ind)*3000; ind=find(mtype>2); mrho(ind)=3300*exp(-2.77e-5*(mTemp(ind)-273)-0.5*0.97e-8*(mTemp(ind).*mTemp(ind)-273*273)-0.32*(1./mTemp(ind)-1/273));
    [gridrho,xgrid,zgrid]=bilingrid(mrho,mx,mz,dx,dz,1,1);
    
    mstressxx=readdata(count,16);
    mstressxz=readdata(count,17);
    mstresstot=(mstressxx.^2+mstressxz.^2).^0.5; clear('mstressxx','mstressxz');
    mstrainxx=readdata(count,18);
    mstrainxz=readdata(count,19);
    mstraintot=(mstrainxx.^2+mstrainxz.^2).^0.5; clear('mstrainxx','mstrainxz');
    mvisc=mstresstot./(2*mstraintot); clear('mstresstot','mstraintot');
    [gridvisc,xgrid,zgrid]=bilingrid(mvisc,mx,mz,dx,dz,1,1); clear('mvisc');

    
    
    ind=find(gridq==0); gridvisc(ind)=nan; gridrho(ind)=nan; gridq(ind)=nan;
    subplot(3,1,1); pcolor(Xgrid,Zgrid,log10(gridvisc)); shading flat; set(gca,'ydir','rev'); axis equal; caxis([17 30]); colorbar; axis([0 2100000 0 710000]); title([num2str(time/(60*60*24*365*1e6)) ' Ma ']);
    subplot(3,1,2); pcolor(Xgrid,Zgrid,gridrho); shading flat; set(gca,'ydir','rev'); axis equal; caxis([2500 3500]); colorbar; axis([0 2100000 0 710000]); title([num2str(time/(60*60*24*365*1e6)) ' Ma ']);
    subplot(3,1,3); pcolor(Xgrid,Zgrid,gridq); shading flat; set(gca,'ydir','rev'); axis equal; caxis([-0.01 0.07]); colorbar; axis([0 2100000 0 710000]); title([num2str(time/(60*60*24*365*1e6)) ' Ma ']);

    print(gcf, '-djpeg100','-r100',['..\outfiles\visc_rho_qz' num2str(count)]);

    
    
 
end
end